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Abstract 

Discrete models of dislocations in cubic crystal lattices having one or two atoms 
per unit cell are proposed. These models have the standard linear anisotropic elas- 
ticity as their continuum limit and their main ingredients are the elastic stiffness 
constants of the material and a dimensionless periodic function that restores the 
translation invariance of the crystal and influences the dislocation size. For these 
models, conservative and damped equations of motion are proposed. In the latter 
case, the entropy production and thermodynamic forces are calculated and fluc- 
tuation terms obeying the fluctuation-dissipation theorem are added. Numerical 
simulations illustrate static perfect screw and 60° dislocations for GaAs and Si. 
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hydrodynamics 
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1 Introduction 



Heteroepitaxial growth is fundamental for manufacturing important nanoelec- 
tronic devices based on self-assembled quantum dots [T1I2] or superlattices [3]. 
Powerful imaging techniques based on scanning probe microscopes have been 
developed which allow to monitor these crystal growth processes down to 
atomic distances. Moreover, these techniques help visualizing defects such as 
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misfit dislocations [45 6 7J (which may act as nucleation sites for quantum 
dots [1]) down to their cores. Understanding epitaxial growth requires under- 
standing and relating processes which cover a wide range of scales [8]. Other 
multiscale problems need to be solved if we try to understand how disloca- 
tions PHU], grain boundaries [3], cracks [11] and other defects control the 
mechanical, optical and electronic properties of the resulting materials |12j . 
Understanding these multiscale processes requires understanding better the 
relation between the defects themselves and observed macroscopic behavior, 
which remains an active area of research. While over mesoscopic length scales 
a coarse grained description based on defect densities makes sense [T3"][T4] , the 
cores of defects need to be resolved on nanometric scales and then coupled to 
the coarser description. Coarse grained descriptions of dislocation dynamics 
include irreversible thermodynamics, as in the works by Holt [13] and by Rick- 
man and Vinals [14], Boltzmann-type kinetic equations, as in the formulations 
by Groma and coworkers [13] (and references therein), or stochastic models 
|16j . The atomic scale near defects can be resolved by ab initio or molecular 
dynamics simulations, which are very costly at the present time. Thus, it is 
interesting to have systematic models of defect motion in crystals that can 
be solved cheaply, are compatible with elasticity and yield useful information 
about the defect cores and their mobility. 

In a previous paper, we have proposed a discrete model of dislocations and 
their motion in cubic crystals with a one atom basis [I7j . In this paper, we 
present an extension of our previous theory to treat crystals with two-atom 
basis in their primitive cells (such as the diamond and zinc-blende structures 
of silicon and gallium arsenide, respectively). Moreover, we explain how to 
include dissipation in the dynamics of the model and how to consider the ef- 
fect of fluctuations by using the ideas of fluctuating hydrodynamics [T8lll9|20j . 
Our model covers length scales in the nanometer range. In principle, to make 
contact with existing mesoscopic theories [T41I15II16] . one should define a dis- 
location density tensor and coarse grain over length scales up to hundreds of 
nanometers. This is outside the scope of the present work. 

The main ingredients entering our discrete model are the elastic stiffness con- 
stants of the material and a dimensionless periodic function that restores the 
translation invariance of the crystal and influences the dislocation size. To be 
precise, consider a simple cubic symmetry with one atom per lattice point. 
Firstly, we discretize space along the primitive vectors defining the unit cell of 
the crystal x = (x, y, z) = (I, m, n)a, in which a is the length of the primitive 
cubic cell, and I, m and n are integer numbers. Secondly, we replace the gra- 
dient of the displacement vector Ui(x,y, z,t) = aui(l,m,n;t) (ui(l,m,n;t) is 
a nondimensional vector) in the strain energy density by an appropriate pe- 
riodic function of the discrete gradient, g(Dj'Ui): We shall define the discrete 
distortion tensor as 
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V1 ? = g{Dfu i ), (I) 
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D 1 m, n;t) — ± [it.j(£ ± 1, m, n; t) — Mj(Z, m, n; £)], (2) 

etc., where g(a;) is a periodic function of period one satisfying g(x) ~ x as x — > 
0. The strain energy density for the discrete model is obtained by substituting 
the strain tensor in the usual strain energy density: 



W = -Cijkieije k i, (3) 
Cii — Ci 

Cijki = C\2 5ij8ki H (8ik$ji + SiiSjk) 

+H (J±JI±JU± _ 5u5 lj 5i k 5 1 i - 5 2i S 2j S 2k 5 2 i - hihjhkhi] , (4) 



fr = 2C44 + C a2 -Cu, (5) 

g(D+Ui) +g(Dfu j 
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(sum over repeated indices is assumed). Here A = C\ 2 , \i = (Cn — C\ 2 )/2 are 
the usual Lame coefficients if H = and therefore the crystal is isotropic. 
Summing over all lattice sites, we obtain the potential energy of the crystal: 



V({ Ui }) = a 3 J2 W(l,m,n;t), (7) 

l,m,n 

in which we have considered the strain energy deensity to be a function of the 
point W(u) = W(l,m,n;t), (l,m,n) = (x,y,z)/a. Next, we find the equa- 
tions of motion with or without dissipation by the usual methods of classical 
mechanics. For conservative dynamics: 
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pa Uiil^m^n^t) 



a dv,i(l, m, n; t) ' 
or, equivalently (see Section l2l) [T7] . 



pa 2 Ui = Yl D j 9'{D+Ui) g(Dfu k )], (9) 

3,k,l 

Here = d 2 Ui/dt 2 and the displacement vector is dimensionless, so that both 
sides of Eq. fl9]) have units of force per unit area. Let us now restore dimensional 
units to Equation (jUJ), so that Ui(x,y,z) = aui(x/a,y/a, z/a), then let a — > 0, 
use Eq. fl9]) and that g(x) ~ x as x — > 0. Then we obtain the usual Cauchy 
equations of linear elasticity: 
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dt 2 



provided the components of the distortion tensor are very small. Far from the 
core of a defect, the discrete gradient approaches the continuous one. Then, 
provided the slope g'(0) is one in the appropriate units, the spatially discrete 
equations of motion become those of the anisotropic elasticity. 

The periodic function g(x) ensures that sliding a plane of atoms an integer 
number of times the lattice distance a parallel to a primitive direction does 
not change the potential energy of the crystal. We choose 



which is periodically extended outside the interval (a — 1/2, a + 1/2) for a 
given a G (0, 1/2). To select a, we calculate the Peierls stress needed to move 
a given dislocation as a function of a and fit it to data from experiments or 
molecular dynamics calculations. Once the discrete model is specified, differ- 
ent dislocation configurations can be selected by requiring that their far field 
should adopt the well-known form of continuous elasticity p2] . 

The rest of the paper is organized as follows. In Section [2J we review the 
derivation of the governing equations with conservative dynamics for simple 
cubic symmetry, and give the numerical constructions of screw and edge dis- 
locations. We use the well known screw and edge dislocations for anisotropic 
elasticity to set up the boundary conditions far from the dislocation core and 
the initial conditions in overdamped equations of motion. Numerical solution 
of these equations yields the static dislocation configuration of our discrete 
elasticity model. In Section [3] we include dissipation and fluctuations in the 
equations of motion. Dissipation is described by a Rayleigh dissipative func- 
tion that is a quadratic functional of the strain rate tensor, which, in turn, 
depends on the discrete distorsion tensor. Since the distortion tensor (con- 
taining finite differences of the displacement vector) and its rate are larger 
near the core of defects, we expect that dissipation will be stronger near the 
core of a moving dislocation than at its far field. Fluctuations are introduced 
via the fluctuation-dissipation theorem and they should be stronger near the 
core of moving dislocations. An extension of our ideas to crystals with more 
complicated symmetries requires formulating our equations in non-orthogonal 
coordinates, which is explained in Section HI The equations of motion for two- 
atom bases are obtained in Section [5] and the corresponding screw and 60° 
perfect dislocations are calculated for diamond and zinc-blende structures. 
Section contains our conclusions. 
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2 Conservative equations of motion for a simple cubic lattice 



In this Section, we shall derive the equations of motion (Q for the conservative 
dynamics given by Firstly, let us notice that 



dW dW de jk 1 d\g(Dfu k ) + g(D+ Uj 

— a jk " 



dui(l,m,n;t) de 3k dui(l,m,n;t) 2 dui(l,m,n;t) 
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5 dui(l,m,n;t) dui(l,m,n;t) 
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where W is a function of the point (/', m', n'), and we have used the definition 
of stress tensor: 

^ = ^' (13) 

and its symmetry, <7jj = cr^. Now, we have 

— — -r[Dju k (l',m',n';t)] = 5 ik {5 IV+1 - 8 a >) 5 mm >5 nn >, (14) 

0Ui{l,m,n;t) 

and similar expressions for j = 2, 3. By using (I12p - dHj), we obtain 



'V ' ' ' ' l,m,n' J 

In this expression, no sum is intended over the subscript i, so that we have 
abandoned the Einstein convention and explicitly included a sum over j. 
Therefore Eq. (jSJ) for conservative dynamics becomes 



pa 2 u t = J2 Dj [<Tij g'(Df Ui )}, (16) 

j 

which yields Eq. (jHJ). Except for the factor g'(D~j~Ui), these equations are dis- 
cretized versions of the usual ones in elasticity [21]. 



2. 1 Static dislocations of the discrete model 



To find the dislocation solutions of our model, we need the stationary solution 
of the anisotropic elasticity equations at zero applied stress corresponding to 
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the same type of dislocation. In all cases, the procedure to obtain numerically 
the dislocation from the discrete model is the same. We first solve the station- 
ary equations of elasticity with appropriate singular source terms to obtain the 
dimensional displacement vector u(x, y, z) = (ui(x, y, z),u 2 (x, y, z),u 3 (x, y, z)) 
of the static dislocation under zero applied stress. This displacement vector 
yields the far field of the corresponding dislocation for the discrete model, 
which is the nondimensional displacement vector: 



U(Z, m, n) = fi((* + Ma,(™ + fl2)a> + W 

(X 

Here 0<5j<l,i = l,2,3, are chosen so that the singularity atx — y — z — 
does not coincide with a lattice point. For a sc crystal, it is often convenient 
to select the center of a unit cell, 5i = 1/2. We use the nondimensional static 
displacement vector U(/,m, n) defined by i fTTj) in the boundary and initial 
conditions for the discrete equations of motion. 

Take for example, a pure screw dislocation along the z axis with Burgers vector 
b = (0,0, b) has a displacement vector u = (0, 0, u 3 (x, y)) with u 3 (x,y) = 
b (27r) _1 tan _1 (?//x) [9j. The discrete equation for the z component of the 
nondimensional displacement u 3 (l,m;t) is: 



pa 2 u 3 = C 44 {£>r [9(Diu 3 ) g'(D+u 3 )) + D 2 \g{D+u 3 ) g'(D+u 3 )}}. (18) 

Numerical solutions of Eq. fTTBl show that a static screw dislocation moves if 
an applied shear stress surpasses the static Peierls stress, \F\ < F cs , but that 
a moving dislocation continues doing so until the applied shear stress falls 
below a lower threshold F c d (dynamic Peierls stress); see Ref. [22] for a similar 
situation for edge dislocations. To find the static solution of this equation 
corresponding to a screw dislocation, we could minimize an energy functional. 
However, it is more efficient to solve the following overdamped equation: 



(3u 3 = C 44 {D-[g(D+u 3 ) g'(D+u 3 )} + D 2 [g(D+u 3 ) g'(D+u 3 )}}. (19) 

The stationary solutions of Eqs. (1T81) and (1T91 are the same, but the solutions 
of f|T9|) relax rapidly to the stationary solutions if we choose appropriately the 
damping coefficient (3. We solve Eq. ffl9l with initial condition u 3 (l,m;0) = 
U 3 (l,m) = b(2?ra)^ 1 tan- 1 [(m + l/2)/(Z + 1/2)] (corresponding to ^ = 1/2), 
and with boundary conditions u 3 (l,m;t) = U 3 (l,m) + Fm at the upper and 
lower boundaries of our lattice. At the lateral boundaries, we use zero-flux 
Neumann boundary conditions. Here F is an applied dimensionless stress with 
\F\ < F cs (the dimensional stress is C44F). For this small stress, the solution 
of Eq. (fl9|) relaxes to a static screw dislocation u 3 (l,m) with the desired far 
field. Figure 3 of Ref. [T7] shows the helical structure adopted by the deformed 
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lattice (l,m,n + u^l^m)) for an asymmetric piecewise linear g(x) as in Eq. 
(jTUl) . The numerical solution shows that moving a dislocation requires that 
we should have g'{D^u^) < (with either j = 1 or 2) at its core [22], which 
is harder to achieve as a decreases. A discusion of the changes in the size of 
the dislocation core and the Peierls stress due to a can be found in Ref. |17j ; 
see in particular Figure 2. Using the same technique, stationary planar edge 
dislocations for an isotropic sc material have been constructed and a variety 
of dipole and loops of edge dislocations have been numerically found [TT] . 



3 Dissipative equations of motion and fluctuations 



3. 1 Equations of motion including dissipation 



Overdamped dynamics obtained by replacing the time differential of the dis- 
placement vector instead of the inertial term in the equation of motion is 
not too realistic. Instead, we can add dissipation to the equations of motion 
by considering a quadratic dissipative function with cubic symmetry: 



R=\C- -r]j y + ve-ik + g ~ ^uSu ~ e^ifak ~ hzhihkf (20) 

For an isotropic body, we have 7 = and then £ and rj are the usual viscosities; 
see Eq. (34.5) in Ref. [21]. The viscous part of the stress tensor is the symmetric 
tensor 



dR _ 

J ik ~. Viklm&lrm V^-U 



ik 



1 / 2 \ rj 
Vmm = 2 l C " 3 V ) 6lk6lm + 2 {6ilSkm + Sim6kl) 

+ 2 { ll km + lm kl _ 5 u 5 lk 5 u 5 lm - 82i5 2 k$2l$2m - 8&5 S k83l83rn) ■ (22) 

In the cubic case, the viscosity tensor rjikim is determined by the three scalar 
quantities (, r\ and 7. For isotropic sc crystals, 7 = 0. Similarly to Eq. ( fTBT) . 
we can show that 



£ R(l\m\n'-t) = -Y.Djl^ l] g\D+u i )], (23) 



dui(l,m,n;t) v ^f n , 





3 



is minus the dissipative force acting on «j. Then the equation of motion in- 
cluding dissipation becomes 
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pa 2 u t = Y, DJ [K- + Eij) g\D+ Ui )} . (24) 



In the isotropic case and taking the continuum limit a — > 0, Eqs. f[24|) with 
fl2TT) and fl22|) yield the viscous Navier's equations for isotropic elasticity 



p w = ^ A " + (A + /i) v( v ■ G) + 71 4? + ( c + 1) v ( v • w) • (25) 



5.5 Entropy production and fluctuations 



Fluctuations may be included in our formulation by using the ideas of Fluc- 
tuating Hydrodynamics [TSPU] . We need to find the entropy production and 
write it as sum of generalized forces and fluxes. Then both the forces and the 
fluxes are identified. The linear relations between forces and fluxes then yield 
the correlations of the fluctuating quantities to be added to the equations of 
motion. 

To find the production of entropy, we need to derive a few formulas. Mul- 
tiplying the conservative equations of motion for the model (Fl6|) by ttj and 
summing, we obtain 



d . ^ 

dt ^ 



^2^u 2 + W(l,m,n;t) 

m,n L i 



(26) 



after some algebra. Provided viscous terms are included in the equation of 
motion, as in Eq. ( 124"1) . we find 



d x ^ 



dt ^ ^2 

l,m,n L i 



E^T^ + W{l,m,n;t) 



E EuiDj^g'iDfui)]. (27) 

l,m,n i,j 



Let Dj4> = 4> — 4>~ i where 4> = <p(l,m,n) is a function of the point (l,m,n). 
(Then <p~ = <fi{l — 1, m, n) for j = 1, and so on). We have ipDjcf) = Dj(ip(j)) — 
<fi~Djtp, and 



E i>dj4>= e djw) - E 4>~ D 74 

l,m,n l,m,n l,m,n 

= E dj{m)- E 



(28) 



Lm,n 



Lm,n 



after a trivial relabelling of indices. Using this formula, Eq. (|2"TI) becomes 
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E 



d_ 

It 



J2^ + W(l,m,n;t) 



i,0 



- E E^W^) (D+ut 

l,m,n i,j 



E E>V<,- (29) 

l,m,n i,j 



where the symmetry of the viscous stress tensor has been used to derive the last 
equality Eq. (1291) describes the production of internal energy due to viscous 
processes. 

If the temperature is not homogeneous, we need to replace the strain energy 
density to leading order by the elastic Helmholtz free energy density: 



F(u; T) = F (T) — (T — T^a^ + - c^e^e^, (30) 

in which the symmetric tensor aij describes anisotropic thermal expansion 
and sum over repeated indices is again implied [21]. Here the material is un- 
deformed at temperature T in the absence of external forces and we assume 
that the temperature change (T — T ) which accompanies thermoelastic de- 
formation is small (linear thermoelasticity) . The stress tensor is now 



&ij = Cijkie k i - aij(T - T ), (31) 

which should be inserted in the equations of motion (Tl6l) or (pMI) . The entropy 
density is S(u; T) = —dF /dT + a^e^ if we ignore the temperature depen- 
dence of the elastic constants. Heat conduction is governed by the equation 
TdS/dt = -dqi/dxi, i.e., 



pc ^t + a - T ^r = qi = (32) 

Here c is the specific heat of the solid and Kij is the symmetric thermal con- 
ductivity tensor. These equations become 



D + T 

pacT + aciijT g'(Djui) Djiii = -D { Q { , Q { = -/fe — — , (33) 

a 

after discretizing. Eq. (|2"9"1) can be rewritten as 



+ W(l,m,n;t) 

Li 



E E ( E ^ y % 

l,m,n j \ i 



DjQ 3 



(34) 
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The right side of this equation is related to the specific entropy (entropy per 
unit mass) s by 



paT 



,ds 
di 



e Dj Qj . 



(35) 



'j 



This can be written as 



D -Q = y s e JL + y qtdt— -y D ;9i 



aT ~3^J ^ V T ^^3 3 aT ^ 3 a j» 

3 *3 3 3 



V 3 1-3 3 

and summing over all points, we find the entropy production: 



E 

Lm,n 



d (P s ) . v D - Qj 
dt ^ 3 aT 



E 

Lm,n 



ij 3 



(36) 



This means that the generalized forces associated to the generalized velocities 
Sy and Qj are —a 3 eij/T and — a 2 Dj~(l/T), respectively. Eqs. £y = r)iji m ei m 
and = (T 2 /a)KijDj~(l/T) then imply that the kinetic coefficients associ- 
ated to Sjj and are k B Tr]iji m / 'a 3 and ks^Kij/a 3 , respectively. Following 
Onsager's ideas as used in Fluctuating Hydrodynamics [T8lll9f20j . we conclude 
that the equations of motion including thermoelastic effects, dissipation and 
zero-mean fluctuations are as follows 



pa 2 Hi 



Dj\(a l3 + E fi + sy) g'{D+ Ui )), 



(sij) = 0, 

(sij(l, m, n; t)s ab (l', mf, n'; t')) = k B T 



T]ijab Tjabij 



(37) 



Su'S mm i5 nn /5(t — t), (38) 



pcf + Ty aijg'iD+m) D+in = — E A r (Q* + 6), 



(6) = o, 

(&(Z, m, n; t)fj(Z', m', n'; t')) = fc B T 



2 n i3 



3 l 



Sii'S mm '5 nn i 5(t — t), 



(39) 



(40) 



with Oij given by fl3T|) . In principle, fluctuations can be included in boundary 
conditions by using the nonequilibrium fluctuating hydrodynamics formalism 
as explained in [23] and in [21] for the case of semiconductor interfaces. In 
crystals with cubic symmetry, the elastic constants and the viscosity tensor 
are given by Eqs. (jlj) and (|22|) . respectively. The thermal conductivity and 
thermal expansion tensors are isotropic, = K8ij, = aSij. Note that the 
correlations of in (|38|) and of in (140]) are proportional to 1/a 3 , which 
becomes <5(x — x') in the continuum limit as a —>■ 0. 
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Note that in our model, dissipation and fluctuations affect all atoms of the 
cubic lattice although we would expect from physical considerations that dis- 
sipation and fluctuations should be more pronounced near the core of moving 
dislocations, as they are directly related to the motion of the atomic con- 
stituents in the core vicinity However our model should also fulfill these ex- 
pectations. Why? Dissipation is described by a Rayleigh dissipative function 
that is a quadratic functional of the strain rate tensor, which, in turn, depends 
on the discrete distorsion tensor. Since the distortion tensor (containing finite 
differences of the displacement vector) and its rate are larger near the core of 
defects, we expect that dissipation will be stronger near the core of a moving 
dislocation than at its far field. Fluctuations are introduced via the fluctuation- 
dissipation theorem and they should also be stronger near the core of moving 
dislocations. 



4 Models of fee and bee crystals with one atom per lattice site 



In this Section, we explain how to extend our discrete models of dislocations 
to fee or bec crystal symmetry, assuming that we have one atom per lattice 
site [17]. For fee or bec crystals, the primitive vectors of the unit cell are not 
orthogonal. To find a discrete model for these crystals, we should start by 
writing the strain energy density in a non-orthogonal vector basis, a±, 02, 03, 
instead of the usual orthonormal vector basis e±, e-i, e?> determined by the cube 
sides. Let Xi denote coordinates in the basis e*, and let x\ denote coordinates 
in the basis eij. Notice that the Xi have dimensions of length while the x\ are 
dimensionless. The matrix T — (ai, a 2 , a 3 ) whose columns are the coordinates 
of the new basis vectors in terms of the old orthonormal basis can be used to 
change coordinates as follows: 

X i = %j X ji X i = ^ij X j- (41) 

Similarly, the displacement vectors in both basis are related by 



u 'i = %j u i = Tiju'ji ( 42 ) 



and partial derivatives obey 



d d d x d 

dx'i 3% dxj' dxi jl dx'j' 

Note that u\ and x\ are nondimensional while Ui and have dimensions of 
length. By using these equations, the strain energy density W = (l/2)cjfc/ m ejfee; m 
can be written as 
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where the new elastic constants are: 



C rspq ~ c ijlmTi r T s j Tl p T qm . (45) 

Notice that the elastic constants have the same dimensions in both the or- 
thogonal and the non-orthogonal basis. To obtain a discrete model, we shall 
consider that the dimensionless displacement vector v! i depends on dimension- 
less coordinates x\ that are integer numbers u[ = u[(l, m, n; t). As in Section 
[2J we replace the distortion tensor (gradient of the displacement vector in the 
non-orthogonal basis) by a periodic function of the corresponding forward dif- 
ference, = g(D^u' i ). As in Eq. (iTTj) . g is a periodic function with g'(0) = 1 
and period 1. The discretized strain energy density is 



W(l,m,n;t) = l -d rspq g{Dtu' r ) g(D+u' p ). (46) 

The elastic constants c' rspq in (]4"5l) can be calculated in terms of the Voigt 
stiffness constants for a cubic crystal, Cn, C44 and C12, which determine the 
tensor of elastic constants (T4]). The elastic energy can be obtained from Eq. 
fj46l) for W by means of Eqs. ([7]). Then the conservative equations of motion 
(ED are 



,<9V ^ lrr x dV 
no = —T l 

p dt 2 n pq du'p 

which, together with Eqs. ([7]) and fj46|) . yield 

<9V- 

P ~Qf = Tiq^pq 1 D j lg'( D t U 'p) C pjrs9(D + U r )} . (47) 

This equation becomes (02) for orthogonal coordinates, = 5{ q /a. 

To add dissipation and fluctuations to these equations, we need to replace 
d pjrs g(D+u' r ) by d pjrs g(D+u' r ) - a' pj {T - T ) + r)' pjrs g'(D+u' r ) D+ii' r + s' pj , in 
which 77' >s is related to the viscosity tensor (1221 in the same way as c' >s 
is related to Ciji m by ([45]) . The random stress tensor s' • has zero mean and 
correlation given by (|38|) with the modified viscosity tensor r/^ ab instead of the 
viscosity tensor (122]) . The heat conduction equations are 
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(^m,n;t)^(l',m\n';t')} = k B T 



pq 



/7--I/7--1 . 



pq 




Su'6 mm f5 nn r6(t-t'), (49) 



(50) 



Note that the both the original and the modified tensors and 
symmetric. 

Once we have derived the equations of motion, stationary dislocations can 
be calculated by first finding the corresponding solution to the equations of 
anisotropic elasticity and using it to set up initial and boundary conditions 
for overdamped equations of motion. For fee and bec crystals, screw and edge 
dislocations have been constructed in Ref. [17] . 



5 Models for diamond and zincblende structures 



Silicon and gallium arsenide are semiconductors of great importance for in- 
dustry that crystalize in the face centered cubic (fee) system. Crystals of these 
materials can be described as a fee Bravais lattice with a basis of two atoms 
per site, which constitute a diamond structure for Si and a zincblende struc- 
ture for GaAs [25] . When growing layers of these materials, defects are very 
important because they act as nucleation sites, and have to be eliminated 
after the growth process has ceased. Among defects, dislocations and misfit 
dislocations are often observed [51T4"] . Thus, it is desirable to have an economic 
description of these defects and their dynamics in terms of control parameters 
such as temperature [12] . A molecular dynamics description is very costly if 
we need to couple atomic details in the nanoscale to a mesoscopic description 
in larger scales that are important in the growth process [26] . In this Section, 
we extend the previous models for cubic crystals with an atom per lattice site 
to crystals having two atoms per site (extension to crystals with more than 
two atoms per site is straightforward). Having two or more atoms per site in- 
troduces new features that are better explained revisiting the classic Born-von 
Karman work on vibrations of a linear diatomic chain[27j. We shall show how 
to obtain the wave equation for acoustic phonons in the elastic limit, directly 
from the equations for the diatomic chain. A similar calculation allows us in 
Subsection 15. 2l to obtain the Cauchy equations for anisotropic elasticity in the 
continuum limit of our discrete models, which are constructed with the aim 
of having exactly this property. Subsection 15.31 shows how to calculate static 
dislocations for GaAs and Si. 
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5. 1 Continuum limit for the linear diatomic chain 



We shall consider a diatomic chain comprising alternatively atoms of masses 
Mi and M 2 whose equilibrium positions are separated a distance a/2. The 
atoms are restricted to move only along the length of the chain. Their dis- 
placement with respect to their equilibrium positions will be denoted by aui 
and avi, respectively, in which / is the cell index, and u\ and vi are dimension- 
less. If (f) is the quadratic potential of interaction between neighboring atoms, 
the equations of motion for the diatomic chain are 



M x ui 



M 2 v x 



vi -ui)a + 



ui - vi-x) a + 



4>" (|) [(vi-ui) - (uj-Ui-i)], 



u l+1 -vi)a + 



yi -ui)a + 



= [{u l+l -v l )-{v l -u l )\. 
If we assume that the solutions of these equations are plane waves, 



(51) 



(52) 



the following dispersion relation is obtained 



(53) 



MiM 2 



[(Mi + M 2 ) =f y (Mi + M 2 ) 2 - 4M X M 2 sin 2 mj\, 



(54) 



in which the minus (resp., plus) sign corresponds to the acoustic (resp., optic) 
branch of the dispersion relation [27]. Moreover, the corresponding amplitude 
ratio for the acoustic branch is 



U 



-M 2 (1 + e - i2 ^) 



V (Mi - M 2 ) - J[M X + M 2 ) 2 - 4MiM 2 sin 2 vrr/ ' 



(55) 



with a similar formula for the optical branch [27] . In the long wavelength limit, 
77 — ^ 0, the acoustic vibrations satisfy 



U — V, u = c 



2-K7] 



(j)"(a/2)a 2 
\ 2 (Mi + M 2 ) 



(56) 
(57) 
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a 2 

In these equations, E and p are the Young modulus and the linear mass density, 
respectively [27J. In the limit as rj — > 0, each cell comprising two atoms moves 
rigidly with a phase velocity c and a wave number 2nr]/a. 

The continuum limit of the diatomic chain equations recovers the acoustic vi- 
brations only. In this limit, I — > oo and a —>■ 0, with fixed x = la. Furthermore, 

a Ui(t) = u(la, t) = u(x, t), a vi(t) = u (la + ^, tj = u (^x + ^, t^j . (59) 

If we now add Eqs. ( I5T!) and (|52|) divide by a, and use (1591) to approximate 
the result, we obtain the following wave equation in the continuum limit: 



d 2 u „ d 2 u , s 

"aF = £ ^- < 60 » 

provided p and are given by Eq. ( |58l) . The wave speed c is then given by 
Eq. ( 1571) . Equation ( l60i) is the elastic continuum limit of the diatomic chain 
equations, which does not contain optical vibrations. 



5.2 Discrete model for a fee lattice with a two-atom basis 



We shall now propose a discrete model for a fee lattice with a basis comprising 
two atoms, of masses Mi and M 2 , respectively. Although this model is much 
more complicated to describe, the key ideas to show that it is compatible with 
anisotropic elasticity are the same as in Subsection 15. II for the diatomic chain. 

The main ideas needed to write a model for this crystal structure are the 
following: 

(1) Write the strain energy corresponding to a fee crystal in a non-orthogonal 
basis with axes given by the usual primitive directions of the fee Bravais 
lattice. 

(2) Write the corresponding strain energy for a fee crystal with two atoms 
per lattice site. 

(3) Restore the periodicity of the crystal by defining the discrete distortion 
tensor as a periodic function (with period 1) of the discrete gradient of 
the displacement vector. 

(4) Define the potential energy of the crystal as the strain energy times the 
volume of the unit cell summed over all lattice sites. Then write down 
the equations of motion for the displacement vectors at each site. 
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(5) Check that the continuum limit of the model yields the usual anisotropic 
elasticity. 



We shall now carry out this program, which is an extension of that presented 
in Section H] for a fee lattice with a single atom per site; see also Ref. [17]. The 
primitive vectors of the fee lattice are 



in terms of the usual orthonormal vector basis ei, 63 determined by the 
cube. From these vectors, we determine the matrix %j to change coordinates 
as in (j4ip and (j42H . In the continuum limit, the strain energy is given by (jUJ) 
with elastic constants given by fj45l) and (T4]). 



Fig. 1. Relevant vectors joining lattice points that are needed to discretize the 
displacement field in a zincblende lattice. All coordinates are expressed in the 
non-orthogonal basis spanned by the primitive vectors a±, a 2 , and 03. (a) The ba- 
sis of a unit cell placed at R' = (l,m,ri) comprises one atom of mass Mi with 
displacement vector and one atom of mass M 2 and displacement vector 

v'^R' + A';i). (b) Discrete gradients involving lattice points closest to R' (resp. 
R' + A') are backward differences from R' + A' (resp., forward differences from R') 
along the primitive directions: Djv'^R' + A';t), (resp., D^u'^R'-jt)), i,j = 1,2,3. 
(c) The auxiliary vectors b\ satisfy a' i + b' i = A', z = 1,2,3. 

Once we have written the strain energy of a fee crystal in the non-orthogonal 
basis spanned by the primitive vectors, we can introduce our discrete model. 
We shall consider a fee lattice with a two-atom basis. In equilibrium, atoms 
with mass M\ will be placed at the lattice sites, so that their displacement 
vectors will depend on integer numbers and time: u' { = «•(/, m, n; t). In equilib- 
rium, atoms with mass M 2 will be placed at the sites of a fee lattice which, in 
the orthonormal basis e^, is rigidly displaced by a vector A = (a/4, a/4, a/4) 
with respect to the first fee lattice; see Fig. [lj In terms of the non-orthogonal 
basis cii, the vector (a/4, a/4, a/4) becomes A' = (1/4, 1/4, 1/4). 




(61) 
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We should define discrete differences of a displacement vector so that a dis- 
crete difference become the corresponding partial derivative in the continuum 
limit. This requirement can be satisfied in more than one way using different 
neighbors of a given lattice point. We shall select only two neighbors of a lat- 
tice point for this purpose, using that the nearest neighbors of an atom with 
mass Mi are atoms with mass M 2 and viceversa. Fig. [1] shows that each atom 
with mass M\ (resp., M2) is linked to its four nearest neighbors having mass 
M 2 (resp., Mi) by A', b\ = A' — a[ (resp., —A', — Thus the nearest neigh- 
bors of an atom with displacement vector v[{Bl + A';t) have displacement 
vectors u'^R'^t) and u'^R' + a'^t), with j = 1,2,3, and the nearest neigh- 
bors of an atom with displacement vector u'^R'^t) have displacement vectors 
v[{R' + A'; t) and v'^R' + A' — a'j] t), with j = 1, 2, 3. These facts motivate our 
definition of discrete differences of a displacement vector. 

Let us define the standard forward and backward difference operators along 
the primitive directions as 



Dff(R') = ±[f(R'±a' j )-f(R')}. (62) 

Then f2^(i?' + A'; t) = D^u'^R'] t) is the discrete gradient of the displacement 
vector v'^R' + A';t) which involves lattice points closest to R' + A', whereas 
Q[j\r'; t) = Djv'^R' +A'; t) is the discrete gradient of the displacement vector 
u'^R'; t) which involves lattice points closest to R'. The distortion tensor of our 
discrete model at the cell R' could be defined as a weighted average of g(Q\f) 

and g(p$), in which g(x) is a period-one periodic function such that g(x) ~ x 
as x — > 0. For simplicity, we shall adopt equal weights in our definition: 



w 



'(#;<) = — J - ' V J (63) 



Obviously, this is reasonable for materials such as Si having a diamond struc- 
ture, and also in the case of atoms of similar size for materials with zincblende 
structure. In the continuum limit a — > 0, the distortion tensor tends to the 
gradient of the displacement vector: 



«,» ~ §| (64) 

In practice, the period-one function g should be fitted to experimental or 
molecular dynamics data, such as the Peierls stress needed for a dislocation to 
move; see Fig. 1 of Ref. [17] for the variation of the Peierls stress as a function 
of the parameter controlling the shape of a piecewise linear function g. The 
positive potential energy of the crystal will therefore be 
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l,m,n 

Here a 3 /4 is the volume spanned by the three primitive vectors. 
The conservative equations of motion are 



(65) 



M 1 



Mo 



dt 2 
d 2 v[ 



tq du' p 

iq pq dv' p 



Using Eq. f[6"5]) . these equations become 



(66) 
(67) 



LUl 0l "' -]%q%~ g lD 7W Pj rs9'{Dfu p ) [g(D:v' r )+g(Dtu' r )}}, (68) 



a" 



dt 2 



4 

4/l/f f) 2 u' 1 

-jT -W = \ T n T M D 1K 3TS g\D- 3 v' v ) [g(Dtu' r )+g(D;v' r )}}. (69) 



If Mi = M.2 (the case of Si), this system of equations is invariant under the 
symmetry: u'^R') <-> —v'^R' + A'). To obtain the continuum limit, we add 
Eqs. fl68|) and fl69|) . take into account the continuum limit (1S"4"|) . and use Eqs. 
fj4TT) to revert to the dimensional orthogonal coordinates. Then the resulting 
equations are those of anisotropic elasticity: 



d 2 Uj d 



dt 2 



dxj 



du r 



(70) 



In this equation, the mass density p is the sum of the masses in the primitive 
cell divided by the volume thereof: 



P 



M 1 + M 2 
a 3 /4 



(71) 



Fluctuations and dissipation can be added to these equations in the same way 
as for the models with one-atom basis. The derivation of these equations in 
non-orthogonal coordinates follows those in Sections [2] and [3l but using the 
following energy instead of ([2] 



d x ^ 

dt ^ 



E \ I ) + E ^ c 'rs Pq g(Dtu r ) g{D+u 



rspq 
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l,m,n i,j 



(72) 



In particular, fl36|) also holds in non-orthogonal coordinates, which then yields 
the same formulas for the fluctuations as in orthogonal coordinates. For a 
zincblende structure, the governing equations including dissipation and fluc- 
tuations are: 



D- U.DfT + S ) , (73) 



a 



2 K ij K ji 



(£(Z, m, n- t)^(l', m', n'- 1')) = k B T 2 -2 £ 5 w 5 mm ,5 nn , 6(t - t% (74) 

J a 

AM, d 2 v' 1 

-jT -W = 2 T * ^ W» + E ^ + 4] (75) 
AM d 2 v f - 1 

at* 1 = 2 T * lT n lD t{Hj + S ^ + 41 ( 76 ) 
(4) = o. 



(4(/,m,n;t)< fe (/',m',n';t')) = 6 w 6 mm ,6 nn , 6(t - 0,(77) 

a pj = C pjab e ab ~ a pj(^ ~ ^o)j ^pj = Vpjab^abi (78) 



together with (jSTjj) . In these equations, = (u^ + wj^)/2 is the symmetric 
part of the distortion tensor ( j63l) . 



5.5 Dislocations in Si and GaAs 



Si and GaAs crystals are face-centred cubic with two atoms per lattice site, 
one at (0,0,0) and the other one at A = (1/4,1/4,1/4). Both atoms are 
identical in Si (diamond structure), whereas they are different in GaAs: one 
atom is gallium and the other arsenic (zincblende structure). Each atom is 
tetrahedrally bonded to four nearest-neighbors, and the shortest lattice vec- 
tor (110)/2 links a second neighbor pair, as shown in Figure [2j The covalent 
bond between two atoms is strongly localized and directional, and this feature 
strongly affects the characteristics of dislocations. In turn, dislocations influ- 
ence both mechanical and electrical properties of these semiconductors [4"][T0"] . 
In this Section, we shall construct straight dislocations for Si and GaAs using 
their elastic stiffnesses and a method for calculating their strain field. 
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Fig. 2. Relevant directions in a fee lattice with a two-atom basis. 

Perfect dislocations have Burgers vectors b = (110)/2 (the same as in the 
case of fee lattices with a one- atom basis) and slip on the close packed {111} 
planes having normal vector n. Their dislocation line vectors £ usually lie 
along (110) directions, forming 60° (perfect 60° dislocations) or 0° (perfect 
screw dislocations) with respect to the Burgers vector. 

We will now construct a perfect screw and a perfect 60° dislocation for GaAs 
lattice, both having Burgers vector b = (1, 0, l)/2. Gliding dislocations having 
this Burgers vector will leave behind a perfect crystal, since it is a primitive 
vector of the lattice [J]. For Si the same construction can be used. 

The tensor of elastic constants is given by (j3J) in terms of the Voigt elastic 
stiffnesses and the degree of anisotropy H of (jSJ). We use the following 
values of the Voigt stiffnesses measured in units of 10 9 Pa at room temperature 
(T = 298 K) [28]: 

C n = 165.6, C12 = 63.98, C 44 = 79.51, for Si; (79) 
C n = 118.8, C12 = 53.8, C 44 = 58.9, for GaAs. (80) 

Between 200 K and 800 K, these constants decrease linearly with increas- 
ing temperature, so that —C^dCij/dT « 10~ 4 K~ , i,j = 1,2,3, [29ll30|31j . 
Such small corrections could be straightforwardly included in our calculations, 
modifying minimally our results. 

To calculate the elastic far field of any stationary straight dislocation, we 
shall follow the method explained in Chapter 13 of Hirth and Lothe's book 
[10] . Firstly, we determine an orthonormal coordinate system e", e 2 ', e'l with 
e'g = — £ parallel to the dislocation line and e 2 ' = n being the unitary vector 
normal to the glide plane. In terms of the new basis, the elastic displacement 
field (u",U2,Us) depends only on x" and on x 2 '. 
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Secondly, we calculate the elastic constants in the reference system e'{, e 2 , e 3 ': 

3 

C ijkl = C ijkl H ~^^{Si n Sj n SknSln $in$ jn$kn$ln) ■ (83-) 
n=l 

Here the rows of the orthogonal matrix S = (e'(, e 2 , e 3 ')* are the coordinates of 
the ef s in the old orthonormal basis ei, en I n the new reference system, 
the Burgers vector has coordinates (b", b 2 , b 3 '). 

Thirdly, the displacement vector {u'{, u 2 , w 3 ) is calculated as follows: 

• Select three roots Pi,P2,P3 with positive imaginary part out of each pair 
of complex conjugate roots of the polynomial det[ajfc(p)] = 0, a^ip) = 

c ilkl + ( c ilk2 + c 'kkl)P + c i2k2P 2 - 

• For each n = 1, 2, 3 find an eigenvector Ak(n) associated to the zero eigen- 
value for the matrix dik{p n )- 

. Solve Re£ti ^(n)L>(n) = bf, and Re^ti Tl =l (c'l 2kl +c'l 2k2 p n )A k (n)D(n) - 
0, in which i = 1, 2, 3, for the imaginary and real parts of D(1),D(2),D(3). 

• For k = 1, 2, 3, u» k = Re[-^ ELi ^(n^n) H< + Pnx' 2 r )}. 

Lastly, we can calculate the displacement vector u' k in the non-orthogonal basis 
di from u k . 

For the perfect 60° dislocation, we have 

e« = -^(1,1,2), <%= +={-1,-1,1), 4= -L(-1,1,0), (82) 

and the resulting dislocation is depicted in Fig. 
For the pure screw dislocation, we have 



e'( = +=(-1,-2,1), e5 = -L(-l,l,l), e - = _L(_i )0 ,-l), (83) 
and the resulting dislocation can be observed in Fig. HI 



6 Conclusions 



We have proposed discrete models describing defects in crystal structures 
whose continuum limit is the standard linear anisotropic elasticity. The main 
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ingredients entering the models are the elastic stiffness constants of the ma- 
terial and a dimensionless periodic function that restores the translation in- 
variance of the crystal (and together with the elastic constants determines 
the dislocation size). For simple cubic crystals, their equations of motion with 
conservative or damped dynamics (including fluctuations according as in Fluc- 
tuating Hydrodynamics) are derived. For fee and bec metals, the primitive 
vectors along which the crystal is translationally invariant are not orthogonal. 
Similar discrete models and equations of motion are found by writing the strain 
energy density and the equations of motion in non-orthogonal coordinates. In 
these later cases, we can determine numerically stationary perfect edge and 
screw dislocations. We have also extended our discrete models to the case of 
fee lattices with a two-atom basis, which includes important applications such 
as Si and GaAs crystals. For GaAs, we have calculated numerically two per- 
fect dislocations which may be used to calculate the structure and motion of 
similar dislocations under stress as explained in Ref. [T7j. Similarly to the case 
of the linear diatomic chain in which there are acoustic and optical branches 
of the dispersion relation, we expect that the dynamics of the discrete models 
with two-atom bases are richer than their continuum limits. 
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Fig. 3. Displacement field in a GaAs lattice created by a perfect 60° dislocation of 
Burgers vector b = (1, 0, l)/2. (a) Reference cubic cell with its eight atoms, (b) One 
layer of a perfect undistorted lattice, (c) The same layer distorted by a perfect 60° 
dislocation. Panels (d), (e) and (f) correspond to (a), (b) and (c), respectively but 
we have depicted only two atoms per reference cubic cell. 



25 




Fig. 4. Displacement field in a GaAs lattice created by a perfect screw dislocation of 
Burgers vector b = (1, 0, l)/2. (a) Reference cubic cell with its eight atoms, (b) One 
layer of cubic cells normal to the Burgers vector for a perfect undistorted lattice, 
(c) The same layer distorted by a perfect screw dislocation. Panels (d), (e) and (f) 
correspond to (a), (b) and (c), respectively but we have depicted only two atoms 
per reference cubic cell. 
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